Quantum critical behavior of the one-dimensional ionic Hubbard model 
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We study the zero-temperature phase diagram of the half-filled one- dimensional ionic Hubbard 
model. This model is governed by the interplay of the on-site Coulomb repulsion and an alternating 
one-particle potential. Various many-body energy gaps, the charge-density-wave and bond-order 
parameters, the electric as well as the bond-order susceptibilities, and the density-density correlation 
function are calculated using the density-matrix renormalization group method. In order to obtain 
a comprehensive picture, we investigate systems with open as well as periodic boundary conditions 
and study the physical properties in different sectors of the phase diagram. A careful finite-size 
scaling analysis leads to results which give strong evidence in favor of a scenario with two quantum 
critical points and an intermediate spontaneously dimerized phase. Our results indicate that the 
phase transitions are continuous. Using a scaling ansatz we are able to read off critical exponents 
at the first critical point. In contrast to a bosonization approach, we do not find Ising critical 
exponents. We show that the low-energy physics of the strong coupling phase can only partly be 
understood in terms of the strong coupling behavior of the ordinary Hubbard model. 
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I. INTRODUCTION 
A. Motivation 

Theoretical studies of the ionic Hubbard model (IHM) 
date back as far as the early seventies (see Ref. [l| and 
references therein). The model consists of the usual 
Hubbard model with on-site Coulomb repulsion U sup- 
plemented by an alternating one-particle potential of 
strength 5. It has been used to study the neutral to 
ionic transition in organic charge-transfer saltsii^ and 
to understand the ferroelectric transition in perovskite 
materialsi^ Based on results obtained from numerical^ 
and approximate methods it was generally believed 
that at temperature T = and for fixed 8 a single phase 
transition can be found if U is varied. This quantum 
phase transition was also interpreted as an insulator- 
insulator transition from a band insulator [U <^ S) to 
a correlated insulator {U ^ 6). In the present paper, we 
discuss in detail how this transition occurs. 

In 1999, Fabrizio, Gogolin, and Nersesyan used 
bosonization to derive a field-theoretical model which 
they argued to be the effective low-energy model of the 
one-dimensional IHMn^ Surprisingly, the authors found, 
using various approximations, that the field-theoretical 
model displays two quantum critical points as U is varied 
for fixed S. For U < Ud the system is a band insulator 
(with finite bosonic spin and charge gaps), as expected 
from general arguments. At the first transition point Ud, 
they found Ising critical behavior as well as metallic be- 
havior in the sense that the gap to the bosonic charge 
modes goes to zero at the critical point only. In the 
intermediate regime, Ud < U < Uc2, a spontaneously 
dimerized insulator phase (in which the bosonic spin and 
charge gaps are finite) with finite bond order (BO) pa- 



rameter was found. The authors argued that the system 
goes over into a correlated insulator phase (in which the 
bosonic charge gap is finite) with vanishing bond order 
and bosonic spin gap at a second critical point Uc2 which 
is of Kosterlitz-Thouless (KT) type. 

Several groups have attempted to verify this phase 
diagram for the IHM using mainly numerical methods. 
Variational and Green's function Quantum Monte Carlo 
(QMC) data obtained for the BO parameter, the electric 
polarization, and the localization length were interpreted 
in favor of a scenario with a single critical point Uc and fi- 
nite BO for U > Uc-^ In a different calculation using aux- 
iliary field QMC, data for the one-particle spectral weight 
were argued to show two critical points with an interme- 
diate metallic phaseiiS Exact diagonalization studies of 
the Berry phaseii and energy gapsiSiiii^ have been in- 
terpreted as favoring one critical poinlii^ or two pointsjii 
in two investigations this issue was left unresolved>i2ii4 
Several density-matrix renormalization group (DMRG) 
studies have been performed focusing on different en- 
ergy gaps, the localization length, the BO parameter, 
the BO correlation function, different distribution func- 
tions, and the optical conduct ivity>i^4i^*iSiii Some of the 
results have been interpreted to be consistent with a two- 
critical-point scenarioii^iiSiii In Ref. Q the signature of 
only one phase transition was found and the possible ex- 
istence of a second transition was left undetermined. The 
phase diagram of the IHM has also been studied using ap- 
proximate methods such as the self-consistent mean-field 
approximation^*'^^'^'', the slave-boson approximationji^ 
and a real space renormalization group method^ Al- 
though these studies led to interesting insights, the va- 
lidity of the approximations in the vicinity of the critical 
region can be questioned on general grounds; therefore, 
we do not focus on these approaches any further here. 



2 



The present situation can be summarized as being highly 
controversial. 

Here we refrain from giving a detailed discussion of the 
merits and shortcomings of the various numerical meth- 
ods used and the possible problems in interpretation of 
numerical results in the literature. Instead, we present a 
detailed study of the T — phase diagram of the one- 
dimensional IHM mainly based on DMRG calculations 
on systems with both open and periodic boundary con- 
ditions (OBC's and PBC's). 

We have calculated a number of different many-body 
energy gaps, including the spin gap, the one-particle gap 
(the energy difference of the ground states with N + 1, 
N, and — 1 electrons), and the gaps to the first ("exci- 
ton" ) and second excited states. A definition of the gaps 
is given in Sec. Ill Al Our results explicitly show that dif- 
ferent gaps associated with charge degrees of freedom do 
not coincide in the thermodynamic limit, although they 
are often believed to in the literature (see also Refs. .16 
and IT^ . Our data show that the exciton gap vanishes 
at a coupling which depends on S and which we define 
as Uci- At this critical point the spin gap remains finite. 
The spin gap vanishes at a second critical coupling, which 
defines our Uc2- 

In addition to the energy gaps, we have determined the 
BO parameter and susceptibility as well as the charge- 
density-wave (CDW) order parameter. Since the single- 
site translational symmetry is explicitly broken due to the 
alternating potential, we will avoid using the term "order 
parameter" in describing the CDW order and instead use 
the term "ionicity" to refer to the difference in occupancy 
between sites on the two sublattices {ua — ns)- We find 
that the ionicity is continuous and non-vanishing for all 
values of the interaction strength. 

From the finite-size scaling of the BO parameter, we 
find a parameter regime with a non-vanishing dimeriza- 
tion starting at Ud and ending at Uc2- We find that the 
transitions at both critical points are continuous. The 
BO susceptibility shows one isolated divergence at Ud 
separated from a region of divergence starting at Uc2- 

We have also investigated the electric susceptibility, 
which is finite in the thermodynamic limit for U < Ud 
and diverges at the lower transition point Ud- For 
U > Ud, the behavior is less clear: there seems to be 
a weak divergence with system size near Uc2 and for 
U > Uc2- This behavior is consistent with that of the 
density-density correlation function, which decays ex- 
ponentially as expected in a band insulator phase for 
U < Ud, but surprisingly decays as a power law with 
an exponent between 3 and 3.5 in the strong coupling 
regime, U > Uc2- 

Using a scaling ansatz for the BO and the electric sus- 
ceptibility we can determine the critical exponents at 
Ud- In contrast to the bosonization approach^, we ob- 
tain critical exponents different from those of the two- 
dimensional Ising model. 

For (almost) all observables, we find that a careful 
finite-size scaling analysis is crucial to obtain reliable re- 



sults in the thermodynamic limit. Furthermore, since it 
is necessary to distinguish between fairly small, but finite, 
gaps and order parameters and vanishing ones, a detailed 
understanding of the accuracy of the DMRG data is es- 
sential. 

In order to obtain a comprehensive picture of the 
ground-state phase diagram, we have studied the differ- 
ent phases (as a function of U) for different (5's which 
cover a wide range of the parameter space. We also con- 
sider the limit of large Coulomb repulsion U oo (for 
fixed 6 and hopping matrix element t) and show that 
some aspects of the physics of the model in this limit 
can be understood in terms of an effective Heisenberg 
model, as has been suggested earlier^ but has recently 
been questioned As a result of our investigations, we 
are able to resolve many of the controversial issues and 
present strong indications in favor of a scenario with two 
quantum critical points. At the appropriate points in the 
paper, we will brieffy comment on the relationship of our 
results with the ones obtained in earlier publications. 

The remainder of the paper is organized as follows. In 
Sec. II Bl we introduce the model and discuss the limits 
in which it can be treated exactly. In Sec. II CI we dis- 
cuss the details of our DMRG procedure. In Sec. m our 
finite-size and extrapolated data for the energy gaps are 
discussed. In Sec. lIIII we present our results for the ionic- 
ity and show that in the large-C/ limit they are consistent 
with analytical results obtained by mapping the IHM to 
an effective Heisenberg model. The BO parameter and 
the related susceptibility are investigated in Sec. IIV Al 
We present results for the electric susceptibility and the 
density-density correlation function in Sec. IIV Bl In the 
numerical calculations of Sees. IIIll to IIVI wc use OBC's, 
for which the DMRG algorithm performs best. To com- 
plete our DMRG study in Sec. El we present results for 
the energy gaps calculated for PBC's and summarize our 
findings in Sec. IVII 

B. Model and exactly solvable limits 

The one-dimensional IHM is given by the Hamiltonian 

where Cj^ {'^ju) destroys (creates) an electron with spin 

a on lattice site j and nj„ — Cj^Cj^. We set the lattice 
constant equal to 1 and denote the number of lattice 
sites by L. Here we study the properties of the half-filled 
system with N — L electrons. 

The system corresponds to the usual Hubbard model 
with an additional local alternating potential. It is useful 
to consider various limiting cases in order to gain insight 
into possible phases and phase transitions. For U = and 
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(5 > 0, the model describes a conventional band insula- 
tor with a band gap 6. Since the alternating one-particle 
potential explicitly breaks the one-site translational sym- 
metry, the ground state has finite ionicity. 

The one-dimensional half-filled Hubbard model with- 
out the alternating potential {5 = 0) and with U > 
describes a correlated insulator with vanishing spin gap 
Ag'^(C/) and critical spin-spin and bond-bond correla- 
tion functionsi2i All gaps associated with the charge de- 
grees of freedom, such as the one-particle gap Af^{U), 
are finitei^^ (The gaps discussed here are defined in Sec. 
Ill Al ) The ionicity and the dimerization are zero for all 
values of U. These two limiting cases suggest that the 
system will be in two qualitatively different phases in the 
limits U ^ 6 and U :$> S. 

In the atomic limit, t — 0, and for < C/ < S, every 
second site of the lattice with on-site energy —6/2 (A 
sites) is occupied by two electrons while the sites with 
energy 6/2 (B sites) are empty. The energy difference 
between the ground state and the highly degenerate first 
excited state is 6 — U. For U > 6, both the A and B 
sites are occupied by one electron and the energy gap is 
U — 6. Thus for t = a single critical point Uc{6) = 5 
with vanishing excitation gap can be found. One expects 
similar critical behavior with at least one critical point 
to persist for the full problem with finite t. 

To describe the physics of the IHM in the limit U » 
t, 6, an effective Heisenberg Hamiltonian 



H 



HB 



J = 



U^-6^ 
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was derived in Ref. Q analogously to the strong-coupling 
perturbation expansion of the usual Hubbard model. It 
has recently been pointed out that this strong-coupling 
mapping does not take into account an explicitly broken 
one-site translational symmetry «ii However, it was shown 
in Ref. Q that the strong-coupling expansion preserves the 
one-site translation symmetry in the effective spin Hamil- 
tonian to all orders in the strong-coupling expansion. In 
addition, the ionicity can be derived directly from the 
effective spin Hamiltonian as follows. The symmetry of 
the Hamiltonian [Eq. (^] implies that after taking the 
thermodynamic limit, nja- — %+2(t for a =t, J, and all j. 
Using the Hamiltonian Eq. and the Hellman-Feynman 
theorem, the ionicity 
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The ground-state energy Eq of the effective Heisenberg 
model [Eq. jSJ] is known analyticalljiSi and, in terms of 
U and (5, is given by 



E^ = L 



ln2 - - 

4 



(5) 



in the thermodynamic limit. In the limit U ^ 6, we can 
thus derive an analytic expression for the ionicity 
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It implies that for any U < oo, the ionicity of the IHM is 
nonzero and for large U vanishes as 1 /U^. Since CDW or- 
der is explicitly favored by the Hamiltonian, it is not sur- 
prising that the ionicity is non- vanishing for all finite U. 
As will be shown in Sec. lIIII this expression shows excel- 
lent agreement with our DMRG data for the IHM. This 
gives us confidence that the effective Heisenberg model 
indeed gives correctly at least certain aspects of the low- 
energy physics. Since the Heisenberg model [Eq. has 
a vanishing spin gapr^ the mapping suggests that the 
spin gap also vanishes in the large-fj limit of the IHM. 

Although the alternating potential breaks the one-site 
translational symmetry explicitly, the model remains in- 
variant to a translation by two lattice sites. This leads 
to a site-inversion symmetry for closed-chain geometries 
with periodic or antiperiodic boundary conditions, a 
symmetry which is not present for OBC's. As pointed 
out in Ref. 0, the ground state of the effective Heisen- 
berg model with periodic boundary conditions for sys- 
tems with An lattice sites or antiperiodic boundary condi- 
tions for systems with An + 2 sites has a parity eigenvalue 
of —1 whereas the ground state for U = has a parity 
eigenvalue of -1-1. This suggests that the IHM under- 
goes at least one phase transition point with increasing 
U for fixed 6. This level crossing will be replaced by level 
repulsion and approximate symmetries for other bound- 
ary conditions. In the thermodynamic limit, the effect 
of the boundaries will disappear and the level repulsion 
becomes vanishingly small. It is important to point out, 
however, that a level crossing on small finite systems does 
not necessarily lead to a first-order transition in the ther- 
modynamic limit; careful finite-size scaling must be car- 
ried out in order to determine the critical behavior. 

From these considerations, one expects to find at least 
one quantum phase transition from a phase with physical 
properties similar to those of a non-interacting band in- 
sulator to a phase with properties similar to those of the 
strong coupling phase of the ordinary Hubbard model. 
However, the details of the transition and the physical 
properties of the different phases remain unclear from 
these arguments. Furthermore, the behavior of the BO 
parameter in the critical region cannot be estimated from 
these simple limiting cases. Therefore, a detailed and 
careful calculation of the characterizing gaps and order 
parameters is necessary. Since no direct analytic ap- 
proach is known to be able to treat the parameter values 
in the critical regime, we restrict ourselves to numerical 
calculations using the DMRG method with the details 
described in the next section. 

In the following, we measure energies in units of the 
hopping matrix element t, i.e., set t = 1. In order to be 
able to cover a significant part of the parameter space, 
we have carried out calculations with ^ = 1, ^ = 4, and 
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6 = 20 for weak interaction values U <^ 5, for strong 
coupling U ^ S and in the intermediate critical regime 
U ^ 6. For the sake of compactness, we will mostly focus 
on (5 = 20 when presenting results that are generic to all 
three J-regimes. 



C. DMRG method 

We have carried out our calculations using the finite- 
system DMRG algorithm. Our investigation focuses on 
the ground-state properties for systems with OBC's, i.e., 
we have performed DMRG runs mostly with OBC's and 
one target state, the case in which the DMRG algorithm 
is most efficient. In order to perform the demanding 
finite-size scaling necessary, we have performed calcula- 
tions for systems with up to L = 768 sites, much larger 
than in an earlier work.^^ 

In order to investigate the low-lying excitations, we 
have also performed calculations targeting up to three 
states simultaneously on systems with OBC's. These nu- 
merically more demanding calculations were carried out 
for systems with up to L = 256 sites for three target 
states and with up to L = 450 sites for two target states. 
In order to compare with exact diagonalization calcula- 
tions and to extend its finite-size scaling to larger sys- 
tems, we have performed calculations for PBC's with up 
to L = 64 sites and one to three target states. In this 
case, the maximum system size is limited by the relatively 
poor convergence. 

The DMRG calculations for OBC's with one target 
state were carried out performing up to six finite-system 
sweeps keeping up to m = 800 states. For more mul- 
tiple states and for PBC's up to 12 sweeps were per- 
formed, keeping up to to = 900 states. In order to test 
the convergence of the DMRG runs, the sum of the dis- 
carded density-matrix eigenvalues and the convergence 
of the ground-state energy were monitored. For OBC's, 
the discarded weight is of order 10^® in the worst case 
and the ground-state energy is converged to an absolute 
error of 10""^ but in most cases the absolute error is 10~^ 
or better. This accuracy in both the energy and the dis- 
carded weight gives us confidence that the wave function 
is also well-converged and that local quantities are quite 
accurate. 

For PBC's, the discarded weight is of the order 10^^ 
in the worst case and the convergence of the ground- 
state energy for most runs is up to an absolute error 
of 10""^ or better, but for extreme cases such as L = 
64 and three target states for parameter values near the 
phase transition points, the convergence in the energy 
is sometimes reduced to an absolute error of only 10^^. 
However, we believe that this accuracy is high enough for 
the purposes of the discussion in section Ivl 

In general, we find that our data are sufficiently accu- 
rate so that extrapolation in the number of states to kept 
in the DMRG procedure does not bring about significant 
improvement in the results (at least for OBC's). Details 



of the extrapolations and error estimates for particular 
calculated quantities are given in the corresponding sec- 
tions. 



II. ENERGY GAPS 

One important way to characterize the different phases 
of the IHM are the energy differences between many- 
body eigenstates. Gaps to excited states can be used to 
characterize phases by making contact with the gaps ob- 
tained in bosonization calculations and also form the ba- 
sis for experimentally measurable excitation gaps, found, 
for example, in inelastic neutron scattering, optical con- 
ductivity, or photoemission experiments. In addition to 
the gaps themselves, however, matrix elements between 
ground and excited states as well as the density of ex- 
cited states are important in forming the full experimen- 
tally relevant dynamical quantities. An example is the 
matrix element of the current operator that comes into 
calculations of the optical conductivity. We have investi- 
gated the behavior of the matrix elements for the dynam- 
ical spin and charge structure factors and for the optical 
conductivity using exact diagonalization on systems with 
both PBC's and OBC's. 

In the following, we present DMRG calculations of the 
gaps to first and second excited states, the spin gap, and 
the one-particle gap in which a careful finite-size scaling 
on systems of up to 512 sites is carried out. As we shall 
see, this is necessary in order to resolve the behavior of 
the gaps in the transition regime and to distinguish be- 
tween scenarios with one or two critical points. 

A. Definition of the gaps 

In this section, we study excitations between a non- 
degenerate 5* = ground state and various excited states. 
In the numerical calculations, we have found that for 
OBC's the ground state is non-degenerate with total spin 
5 = for all parameter values studied here. We define 
the exciton gap 

AE = EiiN,S)-Eo{N,S^O) (7) 

as the gap to the first excited state in the sector with the 
same particle number N and with Sz ^ 0, where Sz is 
the z-component of the total spin. We also calculate the 
expectation value of the total spin operator (S^) so that 
S is known. 

The spin gap is defined as the energy difference be- 
tween the ground state and the lowest lying energy eigen- 
state in the S* = 1 subspace 

As = Eo{N, 5 = 1)- Eo{N, 5 = 0). (8) 

When the first excited state Ei{N, 5) in the Sz = sub- 
space is a spin triplet with 5 = 1, As = Ae. Within the 
DMRG, this gap can be calculated by determining the 
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ground-state energies in different Sz subspaces in two 
different DMRG runs. 

If Ae < As, we call the lowest excitation a charge 
excitation. In fact, exact diagonahzation calculations 
for system with PBC's suggest that the gap Ae corre- 
sponds to the gap in the optical conductivity»i^ We have 
carried out additional exact diagonahzation calculations 
that show that the corresponding matrix elements of the 
current operator are also nonzero for OBC's. We there- 
fore expect that Ae (for excitations with S = and 
when Ae < As) corresponds to the optical gap in the 
thermodynamic limit. ^5, To obtain a deeper understand- 
ing of the excitation spectrum in the critical region, we 
also calculate the gap to the second excited state 

AsE=E2{N,S)-Eo{N,S = 0) (9) 

for selected parameters. 

In the literature, gaps to excitations which can be clas- 
sified as charge excitations are often calculated by tak- 
ing differences between ground-state energies in sectors 
with different numbers of particles (this gap is commonly 
called the "charge gap"). In particular, one can define a 
p-particle gap 

Ap = [Eo{N+p,S:^,J + Eo{N-p,S^,,,J 

~2E„{N,S = 0)]/p (10) 

which is essentially the difference in chemical potential 
for adding and subtracting p particles. The spin 5'^i„ is 
the minimal value, 1/2 or for p odd and even, respec- 
tively. Either the one particle gap Ai or the two particle 
gap A2 are commonly used. The calculation of Ai or 
A2 is numerically less demanding than that of Ae since 
it is sufficient to calculate the ground-state energies in 
the subspaces with the corresponding particle numbers. 
However, since these gaps involve changing the particle 
number and, for p = 1, the spin quantum number, it is 
not a priori clear if they can be used to characterize pos- 
sible phase transition points of the A^-particle system. In 
many cases of interest, the difference between Ai, A2, 
and Ae vanishes for L — > 00, but in other systems (an 
example is the Hubbard chain with an attractive inter- 
action), their behavior differs. As we shall see, Ai and 
Ae do behave differently near Ud- In this work we focus 
our investigation on Ai. We have also calculated A2 and 
find that it behaves similarly to Ai , although it generally 
takes on slightly larger values for finite systems. 

Gaps are also used to characterize the phase diagram 
within the bosonization approach&i^ It is generally be- 
lieved that the bosonic charge gap defined there can be 
identified with the gap to the first excited state with spin 
quantum number 5 = (i.e. the exciton gap Ae [Eq. 0] 
as long as Ae < As) and the bosonic spin gap with As 
[Eq. ©], although a formal proof is missing. 

Based on Ae, As, and Ase and the very limited knowl- 
edge on matrix elements due to the small system sizes 
available to exact diagonahzation, no reliable character- 
ization of the metallic or insulating behavior of different 
phases and transition points can be given. 



B. Gaps to excited states 

In this section, we calculate excited states within the 
Sz — sector. Due to the additional numerical difficulty 
of calculating excited states in the same quantum number 
sector, we are restricted to systems of i = 450 lattice 
sites for Ae and L ~ 256 sites for Ase- 

In Fig. n^a), Ae as a function of U is presented for 
(5 = 20 and various L. For comparison, the spin gap 
for L = 300 is also shown. The exciton gap develops 
a local minimum around U = 21.38, which, for increas- 
ing L, becomes sharper. Furthermore, the value at the 
minimum becomes smaller and seems to approach zero. 
There is a cusp in Ae for all system sizes shown here 
at a certain U to the right of the minimum, and then a 
smooth decay towards zero gap with further increasing 
U. As illustrated for L — 300, this corresponds to a level 
crossing with the first triplet {S = 1) excitation, which 
becomes the first excited state for all larger U values, i.e., 
Ae = As . The data for 6 = 1 and J = 4 behave similarly, 
but the increase to the right of the minimum (up to the 
cusp) is substantially steeper as a function of 6, so that 
it approaches a jump. 

In Fig. ^b) , we display Ae , the gap to the second ex- 
cited state Ase, and the spin gap As (calculated using 
the ground state in the Sz = 1 sector) for L = 128. It 
can be seen that Ase < As for [/-values to the left of the 
minimum in Ae- A similar behavior is found for 6 = 4: 
and i5 = 1. This means that there is more than one S = 
excitation below the lowest lying S* = 1 excitation, con- 
sistent with a scenario in which a continuum of 5 = 
excitations becomes gapless at Ud- This is the scenario 
predicted to occur at the first quantum critical point in 
the bosonization approachi^ Since system sizes for calcu- 
lations of Ase were limited to L = 128 {L = 256 for some 
parameter values), we did not attempt to systematically 
extrapolate Ase to the thermodynamic limit. 

We next discuss the finite-size scaling for Ae to the left 
of the cusp. For U sufficiently far from the critical region 
(i.e., the minimum), the finite-size corrections are small 
and the data can safely be extrapolated to the thermody- 
namic limit using a quadratic polynomial in 1/L, leading 
to a finite exciton gap. Close to the minimum, the scal- 
ing becomes more complicated. At smaller system sizes, 
we find Ae = As and the scaling is nonlinear. How- 
ever, at larger system sizes, there is a crossover to lin- 
ear scaling with Ae{L) ^ As(-/j). The crossover length 
scale becomes larger as U approaches the position of the 
minimum. As a consequence, a reliable finite-size extrap- 
olation in the critical region requires very large system 
sizes. 

To investigate the behavior as L — > 00, we interpolate 
Ae as a function of V for fixed L close to the minimum 
with cubic splines. From the interpolation we can read 
off the minimal value of the gap Amin (i) and the position 
C^min(i) for the different system sizes. Fig. |5| shows the 
resulting Ai„in(i) as a function of XjL for J = 1, 4, and 
20. A linear extrapolation of the data gives Ai„in(L = 
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FIG. 1: (a) The exciton gap Ae for finite system sizes L and 5 = 20. The spin gap As for L = 300 is also shown for comparison, 
(b) The exciton gap, the spin gap As and the gap to the second excited state Ase for L — 128. 



oo,(5 = 1) = 3 X 10"^, A„iin(i = oo,(5 = 4) = 5 X 
10-4, and Amin(i = oo,J = 20) = -1 X 10"''. Within 
the accuracy of our data and our extrapolation, these 
minimal gaps can be considered to be zero. In analogy 
with the atomic limit, we interpret the vanishing of the 
exciton gap as defining a critical point. The critical 
coupling Uci can be determined from fitting UmmiL) to 
a linear function in 1/L, as shown for (5 = 20 in Fig. O 
The extrapolation is similar for the other (5-values and 
we obtain Uci{S = 1) « 2.71, Uci{6 = 4) « 5.6 1, and 
Uci{6 = 20) « 21.39. As wih be discussed in Sec. HVBl 
the vanishing of the exciton gap is accompanied by a 
diverging electric susceptibility. 

In Sec. IIV Al we will present strong evidence in favor 
of a spontaneously dimerized phase for Ud < U < Uc2- 
Since the dimerized phase has an Ising-like symmetry, 
as L oo the ground state in this phase is expected 
to be two-fold degenerate and the exciton gap Ae is ex- 
pected to vanish - at least if the thermodynamic limit is 
taken using PBC's. At first glance this appears to be at 
odds with the increase of Ag as a function of U to the 
right of Uci (but before the cusp is reached) as can be 
observed in Fig. ^a). For finite systems, the OBC's lift 
the degeneracy between the states with the two possi- 
ble bond alternation patterns (strong, weak, strong, . . . 
and weak, strong, weak, . . .), energetically favoring one 
of them which becomes the ground state. We have calcu- 
lated the bond expectation values (see Sec. IIV Ap of the 
ground state and the first excited state on systems of up 
io L — 450 (the largest size we were able to reach) and 
find that the first excited state does not have the oppo- 
site alternation pattern. Instead, the alternation pattern 
is the same as in the ground state near the ends, but 
reverses itself in the middle of the chain. This change in 
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FIG. 2: Finite-size scaling analysis of the minimal value of 
the exciton gap Ae. The solid lines are linear fits through 
the four system sizes shown, L — 256, 300, 350, 400. 



the alternating BO parameter is evenly spread over the 
chain so that it has a cosine-like form with two nodes. 
It is difficult to perform finite-size extrapolation on Ae 
in this region since there are few system sizes and only 
a very limited range of U available. However, one might 
speculate that Ae will remain finite as L — > oo due to 
the pinning of the BO parameter at the ends. 

Sufficiently far from Ud, the data presented in Fig. 
n^a) suggest a linear closing of the exciton gap, which 
gets rounded off in the critical region for finite systems. 
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FIG. 3: Finite-size scaling analysis of the [/-value at the min- 
imum of the exciton gap Ae for 5 — 20. The solid line repre- 
sents a linear least-squares extrapolation of the data yielding 
Uci ^ 21.39. 



The larger L, the closer to Ud the deviation from linear 
behavior sets in. This suggests that Ae ~ Ud — U close 
to but below the first critical point. It implies that the 
product of the critical exponents zivi = 1 at the first crit- 
ical point, where zi is the dynamical critical exponent 
and ui is the exponent associated with the divergence of 
the correlation length. 

Our finding of a vanishing exciton gap at the coupling 
Uci for OBC's is consistent with results obtained using 
PBC's and L = 4n. For this case, a ground-state level 
crossing of two spin singlets at C/ = U^{L, 5) (implying a 
zero exciton gap) was found using exact diagonalization 
of small systems A change of the site inversion 
symmetry at U = Ux was also observed. In Sec. we 
will argue that Ux{L — > (x,S) coincides with Uci{S)}^ 
The presence of the ground-state level crossings might 
lead one to speculate that discontinuous behavior will 
persist in the thermodynamic limit, implying a first order 
phase transition at Ux{L = oo,S). However, we find no 
discontinuous behavior for systems with OBC's, either 
on finite systems or in the L ^ oo extrapolations. In 
order to agree with the results obtained for OBC's in 
the thermodynamic limit, the discontinuous behavior for 
PBC's must become progressively smoothed out as i ^ 
oo. 



C. The spin gap As 



is crucial that the finite-size scaling is carried out care- 
fully and systematically in order to determine the behav- 
ior in the thermodynamic limit. As can be seen in the 
scaling as a function of 1 / L for representative U values in 
Fig.EIb), and as was pointed out in Ref. there is non- 
monotonic behavior as a function of 1/L for U < Ud- In 
addition, the minimum of As as a function of 1/L shifts 
to larger system sizes as the critical region is approached. 
This makes an extrapolation to the thermodynamic limit 
in the critical region a difficult task which requires fairly 
large system sizes. In order to carry out an accurate 
extrapolation, we fit to a cubic polynomial in 1/L. 

Fig-Elshows the extrapolated spin gap for S — 20 pre- 
sented together with the extrapolated values for Ai and 
Ae. All three gaps are approximately equal for U <^ Ud 
(see the inset). Close to the transition, as can be seen 
on the expanded scale in the main plot, Ae goes to zero 
at Ud, while As and Ai stay finite and are (almost; see 
below) equal. For U > Ud, Ae increases until it reaches 
the spin gap As. We find a region of J7 > Ud in which 
As(L = oo) has a value that is clearly nonzero, well above 
the accuracy of the data which is of the order of the sym- 
bol size. The behavior is similar for 6 = 4: (not shown). 
For even smaller values of 6, As close to Ud becomes sig- 
nificantly smaller. As a consequence, the region in which 
As is non-vanishing for U > Ud is less pronounced at 
S = I. In this case, As at Ud is only a factor of six 
larger than the estimated accuracy of our data (this has 
to be compared to the factor of 20 for 5 = 4 and 40 for 
S = 20) with a fast decrease for U > Ud- We take the 
estimate of accuracy from comparison of DMRG calcu- 
lations for the one-particle gap of the usual ID Hubbard 
model with Bethe ansatz results. We find that the dif- 
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worst case. We nevertheless interpret this small spin gap 
to be finite for 3=1 and in a small region oi U > Ud- 
For S substantially smaller than 1, it is impossible to re- 
solve a non-vanishing As at U > Ud using the DMRG. 

The spin gap data in Fig. \S\ indicate that As goes to 
zero very smoothly between 21.55 and 21.8 and remains 
zero from there on. We here define Uc2 as the coupling 
at which As goes to zero. As we have argued in Sec. lIBI 
the mapping onto a Heisenberg model at strong coupling 
[Eq. ©I suggests that the spin gap should vanish at 
sufficiently large U. However, we cannot strictly speaking 
exclude that Uc2 = oo from the spin gap data. We give 
further evidence in support of two transition points at 
finite U below. 

Note that the extrapolated (Fig.[SJ) as well as the large- 
L data (Fig. 0J for As display an inflection point in the 
vicinity of Ud- This might be an indication of a non- 
analyticity related to the phase transition at Ud- 



The spin gap As is shown in Fig. 0] as a function of U 
for S = 20 and system sizes between L = 16 and 512. In 
Fig. ^a), one can see that the spin gap systematically 
scales towards zero above a certain U value. However, it 



D. The one-particle gap Ai 

In Fig. IH^a), Ai as a function of U is shown for S = 20 
and different L. Away from the critical region (which 
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FIG. 4: (a) The spin gap As for finite systems L as a function of U and (b) the finite-size scaling analysis for S — 20 for chosen 
?7-values. The system sizes L — 64, 128, 200, 256, 300, 350, 400, 450, and 512 are shown in (b) and are used for a least-squares 
fit to a third-order polynomial in l/L (solid lines). The dashed line in (b) shows the value of As at the largest system size in 
order to illustrate the non-monotonic behavior. 
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FIG. 5: The exciton gap Ae, the spin gap As, and the one- 
particle gap Ai for 5 — 20 after extrapolating to the thermo- 
dynamic limit L ^ oo. The inset shows the result for a larger 
range of U. 



is between U « 21.15 and U w 22), the finite-L data 
rapidly approach the thermodynamic limit and accurate 
results for Ai{L = oo) can easily be obtained by fitting 
to a polynomial in 1/L. Close to Ud, the data for large L 
develop a minimum. As L increases, the position of the 
minimum shifts to larger {/-values. The shape is quite 
rounded for the small system sizes, but becomes sharper 



for the largest sizes. 

In the critical region, the finite-size scaling is again del- 
icate. We examine Ai as a function of 1/L for a number 
of [/-values near Ud for ^ = 20 in Fig. E^b). The data 
sufficiently away from the minimum (on both sites) shows 
linear behavior in 1/L for smaller system sizes, but then 
deviates from linear behavior and saturates at a finite 
value for larger L-values. This behavior is directly re- 
lated to the L-dependence of the minimum of Ai, which 
shifts to larger U and becomes sharper with increasing 
system size. The scale on which a deviation from the lin- 
ear behavior can be observed shifts to larger system sizes 
as U approaches Ud- In order to perform the finite-size 
scaling analysis, we fit to cubic polynomials in 1/L, as we 
did for the spin gap. We have carried out this procedure 
for (5 = 1 and 4 and find that Ai(L, U) behaves similarly. 

We have extracted the position and value at the min- 
ima by interpolating the data for fixed L with cubic 
splines and then extrapolating to L — > oo with a fit to 
a quadratic polynomial. We obtain f/min('5 = 1) « 2.71, 
Un,in{S = 4) « 5.63, and C/mi„((5 = 20) « 21.40 for the 
positions and Ai{S = 1, C/min) « 0.02, Ai{S = 4, C/min) ~ 
0.05, and Ai{S = 20, C/min) ~ 0.08. The minimal values 
are finite to within the resolution of the data and the 
extrapolation, although the values are small, especially 
at small 6. Therefore, Ai is finite in the critical region 
and is certainly larger than Ae which vanishes at Ud- 
The positions of the minima are very close to, but at a 
slightly larger [/-value than Ud- The largest difference 
C/min(^) - Ud{S) turns out to be 0.02 (for 6 = 4). In Ref. 
[T^ calculation were carried out for 5 = 0.6 (in our units), 
this difference was found to be 0.04, and Ai{5, C/min) was 
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FIG. 6: The one-particle gap Ai for 5 — 20. (a) Results for finite systems with L — 16 through 512. (b) The finite-size scaling 
behavior for L = 64, 128, 200, 256, 300, 350, 400, 450, 512. The solid lines in (b) show least-squares fits to a third-order 
polynomial in 1/L. 



concluded to be zero. The authors interpreted this as an 
indication of a second transition point (in addition to Ud 
which they determined from the vanishing of the exciton 
gap). While we have not carried out calculations at this 
value of S, our results suggest that Ai{5 — 0.6, J/min) is 
(perhaps unresolvably) small, but nonzero. Therefore, we 
believe that C/min is not associated with a second phase 
transition. In fact, as we have seen in Sec. Ill Cl the spin 
gap goes to zero at a substantially higher value of U than 
C/niin, and we associate this value with Uc2- 

Up to a small difference (see Fig. inj Ai(_L = oo) and 
As{L = oo) are equal for U < Ud- In fact, the values 
are virtually identical for the largest few system sizes 
and deviate only at smaller sizes. We therefore believe 
that the difference in the extrapolated gaps stems from 
differences in the fitting to the scaling function at smaller 
system sizes and that Ai(L = oo) = As(i = oo) for 
U < Uynin ~ Uci is cousisteut with our results. At this 
coupling, Ai(L = oo) starts to become larger than Ag 
and as U further increases, grows approximately linearly 
in U as one would expect in a Mott insulator. 

To summarize the behavior of the finite-size extrapo- 
lated gaps, we find that for U <^ Ud, Ae = As = Ai 
as in a non-interacting band insulator. As Ud is ap- 
proached, the gaps to two (or more) S — excitations 
drop below As and at least one of them goes to zero at 
Ud- The one-particle gap Ai reaches a finite minimum 
around Ud and then increases (linearly for large U), and 
the spin gap As goes to zero smoothly at Uc2 > Ud- 
This smooth decay of the spin gap makes it difficult to 
quantitatively estimate Uc2- Since the above behavior is 
similar for the widely different potential strengths stud- 
ied here, (5 = 1, 4, and 20, we believe that it is generic 
for all S. 
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FIG. 7: The ionicity {ua - ub) for 5 = 1,4,20. The solid 
lines indicate analytical results from Eq. @ and the symbols 
numerical DMRG results for L = 32 sites. 



III. IONICITY 

As argued in Sec. II Bl the effective strong-coupling 
model (0) predicts that the ionicity {nA—ns) ^ 1/U^ for 
large U . For t = 0, on the other hand, one expects a dis- 
continuous jump from {ua — ns) = 2 to {ua — ns) = 
at the single transition point Uc- Here we explore the 
behavior of {ua — ng) for all U calculated within the 
DMRG. 
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In Fig. El we compare Eq. © for 5 = 1, 4, 20 and 
various U to results obtained from DMRG with OBC's 
and L = 32. By also considering larger system sizes (up 
to L = 512) and PBC's (up to L = 64), we have verified 
that the L = 32 results shown are already quite close 
to the thermodynamic limit for U 5. On the scale 
of the figure the difference between L = 32 and L = oo 
is negligible. For large U, the DMRG data agree quite 
well with the analytical prediction, Eq. JBJ. This gives a 
strong indication that the large-C/ mapping of the IHM 
onto an effective Heisenberg model^ is applicable at large 
but finite U. It is therefore tempting to conclude that 
Uc2 < oo. One should nevertheless keep in mind that 
the excellent agreement of the numerical data and the 
analytical prediction for the ionicity does not constitute 
a proof of this statement. We will return to this issue. 

The DMRG data for {ua — tlb) for L — 32 shown in 
Fig. 13 are continuous as a function of U for all U. We 
examine (nA — riB) more carefully as a function of system 
size in the vicinity of the first phase transition at Ud for 
^ = 20 in Fig. |H1 The main plot shows DMRG data for 
various L as a function of U for S = 20. While the data 
are continuous as a function of U for all sizes, there is 
significant size dependence between U = 21.2 and 21.5, 
near the first critical point at Ud. We have extrapolated 
the data to the thermodynamic limit using a second order 
polynomial in 1/L and have checked that other extrap- 
olation schemes do not lead to significant differences in 
the extrapolated values. The L = oo extrapolated curve 
is shown in the inset. While the curve is still continuous, 
an inflection point can be observed close to Ud. This 
might be related to non-analytic behavior at Ud- We 
have found similar behavior of {ha — tlb) for 5 — 1 and 
4. 

The behavior of {ua — for PBC's (and L = An), 
which we have checked using the DMRG for up to L = 64, 
is quite different. For finite L, the data display a jump 
discontinuity in the critical region which decreases in size 
for increasing L. The origin of this jump is the ground- 
state level crossing at C/x(L, i5). Since we do not observe 
any discontinuity in the ionicity calculated for OBC's for 
(5 = 1, 4, 20 and up to L = 512, and since the jump 
obtained for PBC's becomes smaller with system size, 
we expect that the jump vanishes in the thermodynamic 
limit and {ua — tib) becomes a continuous function. 



IV. ORDER PARAMETERS AND 
SUSCEPTIBILITIES 



A. The bond order parameter and susceptibility 

The energy gaps have given us indications for two crit- 
ical points. To study the nature of the intervening phase 
and the possibility of dimerization in more detail, we cal- 
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FIG. 8: The ionicity {ua ^ tib) for finite systems with L = 
16, 32, ... , 512 for 5 = 20. The inset shows the L ^ oo 
extrapolated value. 



culate the BO parameter 

(5) = E(-i)^ (4+1.'^.. + ■ (11) 



Since the OBC's break the symmetry between even and 
odd bonds, {B) ^ for all finite systems. Therefore, 
a spontaneous dimerization can be obtained directly by 
extrapolating {B) to L ^ oo, i.e., without adding a 
symmetry-breaking field explicitly. One can form the 
corresponding BO susceptibility xbo by adding a term 



dim 



= pE(-i)' 



c 



(12) 



to the Hamiltonian ^ and taking 



XBO = 



d{B){p) 



dp 



(13) 



p=0 



In practice, the derivative is discretized as [{B){p) — 
{B){—p)]/{2p) where p is taken to be small enough so 
that the system remains in the linear response regime^ 
Due to the additional symmetry breaking by the exter- 
nal dimerization field p, the DMRG runs converge more 
rapidly than in the p — case, making it easier to reach 
larger system sizes. Thus we were able to calculate xbo 
on lattices of up to L = 768 sites. 

Fig. Eta) shows (B) as a function of U for S — 20 
and different L. The data develop a well-defined max- 
imum near Ud for large L. The width of the "peak" 
for L = 512 gives a first indication that there is a re- 
gion in which the dimerization is non- vanishing. Typical 
results for the finite-size scaling of (B) are presented in 
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Fig. Elb). For U ^ Ud, the data extrapolate linearly 
to zero in l/L. In the opposite limit, U ^ Ud, we 
find (B) - with k w 0.5 - 0.6. A similar slow 

decay of the BO parameter has also been found in the 
standard and extended Hubbard models at half-filling.^^ 
The substantial finite-size corrections thus require very 
large systems to distinguish between scaling to zero with 
a slow power-law and scaling to a finite L —>■ oo limit. 
Below, but close to Ud, the data for small L initially dis- 
play power-law- like finite-size scaling with k < 1, but for 
larger system size, one finds a crossover to a linear scaling 
of the BO parameter (to zero) as L — > oo. There is also 
a crossover in the behavior for t/- values near but above 
Ud- One again finds a crossover from a power law with 
K < 1 for smaller system sizes to linear behavior that 
can be extrapolated to finite values of (i?)oo for larger 
system sizes. The crossover length scale increases as U 
approaches Ud until it becomes larger than the largest 
system size considered here. This length scale can be 
used to estimate the correlation length, which diverges at 
the first (continuous) critical point. We have been able 
to calculate Lc for U -values on both sides of Ud and find 
that it diverges approximately linearly in | [/ — Ud \ ■ This 
implies vi = 1 (see also below). Taking into account that 
zivi ~ 1 as extracted from the linear closing of Ae, one 
finds zi — 1 for the dynamical critical exponent. 

This diverging crossover length scale makes it essential 
to treat system sizes that are significantly larger than the 
scale Lc, even close to the critical point Ud- In order 
to obtain reliable results, we have calculated {B) l for a 
number of system sizes L > 200. In carrying out the 
finite-size extrapolation, we fit to a linear form for the 
largest system sizes if it is clear that has been reached, 
as can be seen in the inset of Fig.|5Jb). 

In Fig. ^1 the finite-size extrapolation (B)^ is shown 
as a function of U for 5 = 1, 4, and 20. As can be seen, 
{B}^ = to well within the error of the extrapolation 
for U < Ud- For U > Ud, we find a region of width be- 
tween 0.2 and 0.4 (i.e., a factor of 5 to 10 larger than the 
extent of the dimerized phase claimed to be found in Ref. 
IT^ in U in which (B)^ is distinctly finite. The onset of 
finite (B)^ at Ud is rather steep for all three values of 6, 
but seems to be continuous. This steep onset suggests a 
critical exponent of the order parameter that is substan- 
tially smaller than 1. Within bosonization the first crit- 
ical point was predicted to be Ising-like with /3i = l/Sn^ 
The fall-off to zero as U increases, on the other hand, 
is slow, with a small or vanishing slope. This behavior 
would be consistent with a second critical point at which 
the critical exponent for the order parameter is larger 
than one or at which a higher order phase transition such 
as a Kosterlitz-Thouless transition takes place. ^ As can 
be seen by comparing Fig. llOf a). (b), and (c), the height 
of the maximum increases with increasing S- For S signif- 
icantly smaller than 1, the BO parameter is so small that 
it cannot be concluded to be finite within the numerical 
accuracy of the DMRG. For the couplings at which the 
finite dimerization sets in we obtain Ud{S = 20) w 21.39 



and Ud{S = 4) « 5.61, which are in excellent agreement 
with the results obtained from the vanishing of Ae. The 
value obtained for 6 = 1, Ud{6 = 1) ~ 2.67, is also 
in reasonably good agreement with the results obtained 
from the analysis of the gaps. 

While our data suggest that a critical coupling Uc2, 
with (B)^ = ioi U > Uc2, exists, no reliable quanti- 
tative estimate of Uc2 can be given based on the DMRG 
data for the BO parameter. Due the close proximity of 
the two critical points, we were not able to obtain quan- 
titative results for the critical exponents /3i and f32 at the 
critical points, either by a direct fit of the L = oo results 
or by a scaling plot of the finite-size data. As discussed 
next, accurate exponents at Ud can be extracted from 
both the BO and the electric susceptibilities, and a more 
accurate estimate of Uc2 can be obtained from the BO 
susceptibility. 

In order to understand the behavior of the BO suscep- 
tibility, it is useful to first examine the behavior of the 
BO parameter (B) as a function of the applied dimer- 
ization field p- This quantity is shown in Fig. 1111 for 
6 = 20, three representative values of U, and different 
system sizes. For U = 19 < Ud, the system is in a phase 
with vanishing BO parameter, and the slope at p = re- 
mains finite for all system sizes, corresponding to a finite 
susceptibility. The value U = 21.42 is in the intermedi- 
ate regime where we have found a finite BO parameter 
in the thermodynamic limit. As can be seen in the main 
part of the figure, a jump in {B){p) develops. As the 
system size increases, the absolute value of dimerization 
field at which the jump occurs becomes smaller. This 
is the behavior expected in a dimerized phase in a sys- 
tem with OBC's. Therefore, the jump in {B){p) provides 
additional evidence in support of an intermediate phase 
with finite dimerization. For the approximate calculation 
of the susceptibility xbo « [{B){p) - {B){-p)]/{2p), we 
have taken p = 10^^ which is small enough to stay to the 
right of the jump for all system sizes considered. Finally, 
for [/ = 50 ^ Uc2, {B){p) goes to zero for \p\ and 
increasing system size indicating a phase without spon- 
taneous dimerization. However, the slope at small \p\ 
becomes steeper with increasing system size, indicating 
a divergence of xbo ■ 

In Fig. El the BO susceptibility as a function of U 
is shown for (5 = 1, 4, 20 and different L. For all 5 val- 
ues, one observes a two-peak structure that becomes pro- 
gressively more well-defined with increasing system size. 
There is a narrow peak at a [/ -value that agrees well with 
Ud determined earlier whose height grows rapidly with 
system size. It signals the onset of spontaneous dimeriza- 
tion. For somewhat larger U there is a minimum in xbo , 
surrounded by narrow region in which its value seems to 
saturate with system size. For still larger [/-values, a 
second, broad peak develops. The position of this sec- 
ond maximum is roughly at Uc2, the [/-value at which 
the BO parameter vanishes. We argue that the second 
peak is related to the second phase transition from the 
dimerized phase into an undimerized phase. To the right 
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FIG. 9: (a) The bond order parameter {B) for 5 = 20 for finite systems as a function of U for various system sizes, (b) The 
scaling of the data as a function of the inverse system size 1/L. The solid lines are least-squares fits to the data as described 
in the text. The inset shows an expanded view of the scaling for U values near the critical point. 



of the second peak, xbo does not seem to saturate for in- 
creasing system size, implying that Xbo is divergent for 
all U > Uc2- One can understand this divergent behavior 
by studying the BO susceptibility for the ordinary Hub- 
bard model Xbo • One finds that Xbo divergent for 
all J7 > because the bond-bond correlation function is 
criticalr^i A finite-size extrapolation of xbo is shown in 
Fig- El for large [/-values for both (5 = and 5 — 20. We 
find a power law divergence, xbo(-^) L'^ , with C, ~ 0.68 
for the ordinary Hubbard model and C, ~ 0.65 for the 
IHM. These values are in good agreement, considering 
the accuracy of the fit and additional finite-size effects. 

Since xbo diverges for all U to the right of the sec- 
ond peak, it is difficult to accurately determine the crit- 
ical coupling Uc2- However, two different ways of esti- 
mating 17^2(5) under- and overestimate its value. In the 
first method, Uc2 is estimated as the lowest U-value for 
which Xbo seems to diverge for increasing L and the 
available system sizes. It is then still possible that there 
is a crossover above a length scale unreachable by us and 
Xbo scales to a finite value. 

This tends to underestimate Uc2- In the second 
method, Uc2 is taken to be the position of the second 
peak at fixed L, extrapolated to L ^ 00. Since the peak 
position decreases for increasing L, this method tends 
to overestimate Uc2- From these two procedures, we ob- 
tain the bounds 21.55 < Uc2{S = 20) < 21.69. For the 
other values of 6, it is very difficult to accurately de- 
termine the lower bound with the data available. We 
therefore only give the upper bound Uc2{S = 1) < 2.95 
and Uc2{6 = 4) < 5.86. 

It is generally believed that a quantum critical point is 
accompanied by a vanishing characteristic energy scale. ■^^ 



At Uc2 the most obvious candidate is As , consistent with 
our numerical data (see Figs. 0] and 0) and implying that 
Uc2 = Uc2- This is assumed in the following discussion. 

Since the peak in xbo at Ud is well-defined and has a 
clear growth with system size, it is reasonable to perform 
a finite-size scaling analysis. We use a scaling ansatz of 
the form 

x(C/,i)=i2-''x(i/0, (14) 

with ^ - |C/ - Ucl'". As can be seen in Fig. [T^dl. 
data for S = 20 and system sizes of L = 128 and greater 
collapse onto one curve. The best fit is obtained with 
Uci = 21.385 and the critical exponents 771 = 0.45 and 
vi — 1. The latter value is consistent with the value 
i>i = I extracted from the divergence of the length scale 
discussed above. Note that the value of 771 is not in agree- 
ment with the value expected in the two-dimensional 
Ising transition, ry = l/4i^ We have also applied the scal- 
ing ansatz for S — I and 4. For decreasing 6, the quality 
of the collapse of the data for the available systems sizes 
becomes poorer and the extracted exponents therefore 
become less reliable. The best fit is again obtained with 
i^i = 1 for both S, 7n{S = 4) 0.55 and r/i((5 = 1) « 0.65. 
For the critical coupling we obtain Uci{5 = 1) « 2.7 and 
Uci{S = 4) w 5.6, in excellent agreement with the values 
found by other means. 

It is also possible to collapse the finite-size data onto 
one curve at the second transition point using the scaling 
ansatz (I14f) . We find that the best results are obtained for 

~ exp{A/{U — U'^2)^)i indicating that the divergence 
of the susceptibility at Uc2 may indeed be exponential 
as expected for a KT-like transition. However, fitting 
the limited amount of data available to this form does 
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FIG. 10: The bond-order parameter {B)aa in the thermody- 
namic limit for (a) 5 = 1, (b) 5 = 4, and (c) 5 = 20 plotted 
as a function of U near the transition points. 



not produce completely unambiguous results for all fit 
parameters. 

Therefore, we have not further attempted to obtain 
results for A, B, [/^2i ^^^d 772 with this method. 



B. The electric susceptibility and the 
density-density correlation function 

In order to further investigate the physical properties 
of the different phases and transition points, we calculate 
the electric polarization and susccptibilityn^ The polar- 
ization is given by 



{P) 



-E 

3 



Xj {riji^+nji) , 



(15) 



where Xj = j — i/2 — 1/2 is the position along the chain, 
measured from the center. The polarization is the re- 
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FIG. 11: The BO parameter {B) as a function of applied 
dimerization field p for 5 = 20 and U = 21.42. The upper 
inset shows data for [/ = 19 and the lower inset data for 
U = 50. 



ponse due to a linear electrostatic potential 



(16) 



which is added to the Hamiltonian The electric sus- 
ceptibility 



Xcl 



d{P){E) 



dE 



(17) 



E=0 



is the susceptibility associated with this field. 

The electric susceptibility has been used to investi- 
gate the metal-insulator transition in the t-i'-Hubbard 
modeli^ In this model, both a phase in which Xc\ diverges 
as (a perfect metal) and a phase in which for increas- 
ing system size Xc\ scales to a finite value (an insulator) 
were found when varying U for fixed nearest-neighbor 
hopping t and next-nearest-neighbor hopping t' . 

In contrast to the ordinary Hubbard model, the po- 
larization does not always vanish at field i? = in the 
IHM. For = 0, (5 > 0, one finds (P) = -1/2. This 
is due to the alternating ionic potential which induces a 
charge displacement to the sites with lower potential en- 
ergy. Due to the OBC's, a chain with even length L starts 
and ends with a different potential, inducing a dipole mo- 
ment. This is a boundary effect. In the strong coupling 
limit, [/ ^ (5, we find that (P) — > 0, as expected. The 
electric suceptibility Xci can be calculated by discretizing 
the derivative as [(P)(P) - (P)(P = 0)]/P. The field 
E must be taken to be small enough so that the system 
remains in the linear response regime»2£ Note that it is 
necessary to subtract (P)(P = 0) since it is nonzero in 
general. 
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FIG. 12: The BO susceptibility xbo as a function of U for (a) 5 — 1, (b) S — 4, and (c) 5 = 20 and different L. (d) A scaling 
analysis of the S — 20 data from (c). 



A plot of Xci as a function of U for various system sizes 
is shown in Fig. I14r al for S — 20. For U <^ Uci and in- 
creasing L, Xc\ converges to a finite value, similar to the 
behavior in a non-interacting band insulator and in the 
correlated insulator phase of the t-t'-Hubbard modeli^S 
The data clearly develop a maximum at Uci whose height 
increases markedly with system size, indicating a diver- 
gence at the first critical point. The finite-size scaling 
of this height is consistent with a power-law increase, 
i^-j)!^ with rji « 0.46. This increase is weaker than 
the divergence (which implies r/ = 0) found in Ref. 
Isol and associated with a perfect metal. For U slightly 
larger than Ucl^ the data again seem to saturate with 
system size. Assuming the scaling form of Eq. (I14|l . the 
data close to Uci can be collapsed on a single curve as 



demonstrated in Fig. I14r b). The best fit is obtained for 
vi = 1 and rji « 0.45. Both of these exponents are 
in excellent agreement with those found in the scaling 
analysis for xbo- We have carried out a finite-size scal- 
ing analysis for 5 = 4 and 6 = \ and also find diverg- 
ing peaks at ?7ci, as well as collapse of the data onto a 
single curve using the scaling form 114(1 with exponents 
rii{5 = 1) = 0.52, rii{5 = 4) = 0.45, and = 1 (for 
both 5). The critical [/-values obtained from this scaling 
procedure are Uci{5 = 1) = 2.68, Uci{5 = 4) = 5.59, and 
Uci{5 = 20) = 21.38, which compare well to the values 
for the critical coupling obtained from the gaps and from 
the BO parameter and susceptibility. 

The data for 5 = 20 and J = 4 for the largest sys- 
tem sizes, L = 256 and L = 512, suggest that a second 
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FIG. 13: The BO susceptibility xbo as a function of 1/L 
for the ordinary Hubbard model {5 — 0) and U — 10 and the 
ionic Hubbard model for 5 = 20 and U = 50. DMRG data are 
indicated by the corresponding symbols and the solid curves 
represent a least-squares fit to the indicated forms. 

peak may develop around Uc2- In order to investigate 
the behavior of Xci{L) more precisely in this region, we 
fit a quadratic polynomial to (P) (E) through several data 
points and then take the derivative of this fit function at 
E — 0. This procedure should eliminate errors caused 
by a small linear response regime. Results obtained from 
this procedure for i5 = 20 indicate a weak divergence at 
U = 21.65, corresponding to a [/-value near Uc2- In addi- 
tion, we find an even weaker divergence for all U > 21.65. 
The larger the U -value, the smaller the coefficient of the 
diverging part, so that the divergence is very difficult to 
observe numerically deep in the strong-coupling-phase. 
One generally expects the divergence of Xci to be con- 
nected to the closing of a gap to excited states which 
possess at least some "charge character" (in the sense 
discussed below). At Ud the divergence is accompanied 
by the closing of the exciton gap, leading to a consistent 
picture. 

The situation is less clear for U > Uc2- This issue can 
further be investigated by examining the behavior of the 
density-density correlation function 

Cdcn{r) = {n^rii+r) - (rii) (rii+r) , (18) 

shown in Fig. E| for 5 = 1 and different U > Uc2- Here 
we have averaged over a number of i-values (typically 
six) for each r and have performed the calculation on an 
L = 256 lattice. For each value of U, it is evident that 
the correlation function behaves linearly on the log-log 
scale above some value of r, indicating that the domi- 
nant long-distance behavior is a power law. (For r close 
to the system size, finite-size effects from the open bound- 
aries also appear.) Note that the sign of the correlation 



function is negative for r > 0, so that the negative is 
plotted. A least-squares fit to the linear portion of the 
curve yields an exponent of approximately 3 — 3.5 for all 
values of [/ > Uc2- This behavior is markedly different 
from the behavior for U < Ud, where we find a clear 
exponential decay as in a non-interacting band insulator, 
and from the behavior at Ud, where we find a power law 
decay with an exponent of « 2. Note that if the decay 
were exponential for U > Uc2, we would expect the cor- 
relation length to change quickly with U , leading to a 
marked variation in the slope. We have ruled out finite- 
size effects as an origin of the power-law tails as well as 
possible symmetry breaking due to the OBC by compar- 
ing calculations for L = 128 and L = 256 with OBC and 
L — 64 with PBC, which yield identical values except for 
distances r near the lattice size (or half the lattice size 
for PBC's). 

We have performed calculations for i5 = 20 and find 
similar behavior. The exponent of the power-law tails has 
a comparable value to the ones given above, even at very 
large [/-values such asU = 50, where the prefactor of the 
power-law part is « 2 x 10~^. It therefore seems justified 
to conclude that this power-law decay is a generic feature 
of the strong-coupling phase for all 5. 

Our findings for Xci and CdGn(?') are consistent with 
a scenario in which there is a continuum of gapless ex- 
citations for U > Uc2, where matrix elements of charge 
operators such as the density rij = riji + riji, are non- 
vanishing for some of the states belonging to this con- 
tinuum. These are the states mentioned above which 
possess charge character. To further confirm this idea, 
we have calculated matrix elements (m| Uj |0), where |m) 
denotes the m-th excited state and |0) the ground state, 
for up to m = 4, (5 = 20, U > Uc2, and L = 32. We 
find that the third excited state is the first S ~ state, 
both for the ordinary Hubbard model and the IHM (the 
fourth state as well as the m = 1,2 states have S = 1). 
For the ordinary Hubbard model, {3\nj |0) vanishes for 
all j to within the accuracy of our data and this S — 
state can be classified as a spin excited state since its 
excitation energy is well below the charge gap. In con- 
trast, (3|nj |0) is nonvanishing for the IHM and shows 
a non-trivial dependence on j which has a wavelength 
of approximately the lattice size, implying that the wave 
vector characterizing the excitation is near q = Q. 

As a consequence, this state contributes to the dy- 
namical charge structure factor in the IHM but not in 
the ordinary Hubbard model. This shows that although 
several similarities between the strong coupling phase of 
the IHM and the Hubbard model were found, low-lying 
excitations in both models are of quite different nature. 
As we have verified, the energy of |3) becomes smaller 
for increasing U, in contrast to the behavior of the one- 
particle gap which increases linearly with U. Due to the 
numerical effort necessary to target such a large number 
of states, we were unable to perform these calculations 
on larger lattices in order to carry out a finite-size scaling 
analysis of the matrix elements. 
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FIG. 14: (a) The electric susceptibility Xci for S = 20 plotted as a function of U. (b) A scaling analysis of the data of (a) 



It is important to note that the power-law decay of 
Cdon('') and the divergence of Xei for U > Uc2 do not 
necessarily imply that the Drude weight is finite in this 
parameter regime or near Ud, where Xoi diverges roughly 
as L^'^ . Therefore, we refrain here from classifying the 
dimerized phase, the U > Uc2 strong coupling phase, and 
the transition point Ud as being metallic or insulating. 
For U < Uci all our results are similar to those found in 
a non-interacting band insulator. To further investigate 
the metallic and insulating behavior in different parts of 
the phase diagram, it would be necessary to calculate dy- 
namical correlation functions using, e.g., the dynamical 
DMRG. Such an investigation would exceed the scope of 
the current paper. 



PERIODIC BOUNDARY CONDITIONS 
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Up to this point, we have only presented DMRG results 
obtained for systems with OBC's. Here we will argue 
that the results for energy gaps determined for PBC's are 
consistent with the ones discussed in Sec.|n] We present 
further evidence that Uc2 obtained from the closing of 
the spin gap for OBC's and the coupling constant Uc2 at 
which the BO parameter vanishes coincide. 

Here we investigate the level crossing point of the 
ground state and the first excited state U^, as well as the 
crossing of the first and second excited states C/xx- In Ref. 
[Til (and further references therein) , the crossing points of 
ground and excited states of finite systems are associated 
with phase transition points. In particular, the crossing 
of the ground state and the first excited state with op- 
posite site-inversion parity were shown to correspond to 
a jump in the charge Berry phase. The crossing of the 



FIG. 15: The negative of the density-density correlation func- 
tion — [{niUi+r) — (rii) (rii+r)] for U > Uc2- The indicated 
lines are least-squares fits over a range of r in which the be- 
havior is linear on the log-log scale. 



first and second excited states which are spin singlets and 
spin triplets with opposite site-inversion symmetry and 
zero total momentum was argued to be associated with 
a jump in the spin Berry phase and to characterize a sec- 
ond transition. While the direct calculation of the Berry 
phases is beyond the scope of this paper, it is possible for 
us to analyze the finite-size behavior of the level cross- 
ings C/x and C/xx- We therefore must calculate the ener- 
gies of the ground state and the first two excited states 
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FIG. 16; Crossing points of excited states for PBC's for 5 — 1 
as a function of inverse system size. The inset shows the 
extrapolation of As(f/xx) as a function of inverse system size. 



simultaneously. Sufficiently accurate DMRG results for 
these energies can only be obtained for system sizes of 
up to L = 64, small compared to the ones studied for 
OBC's, but nevertheless much larger than the ones that 
can be reached with exact diagonalizationiiiiiSiiiii We 
show results for PBC's for systems with L = 12, 16, 24, 
32, and 64, i.e., system sizes with An sites, so that the 
site-inversion symmetry of the ground state is guaranteed 
to change sign with U, as discussed in Sec. II Bl We find 
non-monotonic behavior of the level crossing points as a 
function of system size at a scale beyond the system sizes 
which can be investigated using exact diagonalization. 

The finite-size scaling of U-^ and t/xx for (5 = 1 is shown 
in Fig. 1161 The error bars result from the uncertainty 
in determining the closing or crossing points as well as 
from the poorer convergence of the DMRG algorithm for 
PBC's. Nevertheless, ioi 6 — I they only arc of the order 
of the symbol size or smaller. For 5 = 4 and S = 20, the 
convergence for L = 64 is poorer around and U^x, but 
we obtain qualitatively similar behavior up to the larger 
error bars. 

For all S studied, the finite-size extrapolation of 
leads to critical couplings in agreement with the ones 
given in Sec. Ill CI up to the smaller numerical accuracy 
available with PBC's. Using a quadratic polynomial for 
the extrapolation, we find Ux{S = 1) = 2.71, Ux{S = 4) = 
5.63 and Ux{S = 20) = 21.42. (Due to complicated finite- 
size effects, we only use the data for L < 32 for 5 = 4 and 
6 = 20.) The angle of crossing of Eq and Ei decreases 
with increasing system size. This is consistent with a 
continuous critical behavior at Ud in the thermodynamic 
limit. 

The non-monotonic behavior of U^x, as seen in Fig. 
1161 makes an i ^ oo extrapolation difficult. In fact, 
an extrapolation using the system sizes available to ex- 
act diagonalization would give a Uc2 which is substan- 



tially larger than if the two largest system sizes were in- 
cluded. This could explain the discrepancy in the size 
of the region between the two critical points found here 
and obtained in Ref. By extrapolating the finite- 
size data using a quadratic polynomial in 1/L, wc obtain 
;7xx(i = Qo,S ^ 1) ^ 2.84, Uxx{L = oo,S ^ A) ^ 5.97, 
and Uxx{L = oo,S = 20) = 21.75. The values for Uxx and 
for Uc2 obtained from the BO susceptibility are in fairly 
good agreement. The differences indicate that even larger 
system sizes are needed to perform an accurate finite-size 
extrapolation for PBC's. 

Since the transition Uc2 is associated with the clos- 
ing of the spin gap, the gap to the excitations at Uxx 
should scale to zero with system size. The inset of Fig. 
1161 shows that Ae(L'xx) = Ase(Kcx) indeed closes in the 
thermodynamic limit. Since one of the two states that 
are degenerate at Uxx is a spin triplet, this implies a 
vanishing spin gap. We thus obtain further (indirect) ev- 
idence that the couplings at which the BO susceptibility 
diverges and the spin gap closes coincide, consistent with 
a two-critical-point scenario. In particular, the angle of 
the crossing of the first and second excited state also de- 
creases with increasing L, consistent with a continuous 
transition at Uc2- 



VI. SUMMARY 

In this paper, we have presented density-matrix renor- 
malization group results that elucidate the nature of the 
quantum critical behavior found in the half-filled ionic 
Hubbard model. By carrying out extensive and pre- 
cise numerical calculations and by carefully choosing the 
quantities used to probe the behavior, we have been able 
to investigate the structure of the transition more accu- 
rately than in previous work. This has allowed us to 
resolve a number of outstanding uncertainties and ambi- 
guities. We have worked at three different strengths of 
the alternating potential S covering a significant part of 
the parameter space and find the same qualitative behav- 
ior for all three J-values. In particular, we have carried 
out extensive finite-size scaling analyses of three differ- 
ent kinds of gaps: the exciton gap, the spin gap, and 
the one-particle gap. We find that for fixed S and in the 
thermodynamic limit, the exciton gap goes to zero as a 
function of J7 at a first critical point Ud, the spin gap 
goes to zero at a distinct second critical point Uc2 > Ud 
and is clearly nonzero at Ud- The one-particle gap (the 
two-particle gap behaves similarly) reaches a minimum 
close to Ud, but never goes to zero and never becomes 
smaller than the spin gap. 

Due to the explicitly broken one-site translational sym- 
metry, the ionicity is finite for all finite U. For U ^ 6 the 
ionicity found numerically agrees very well with the one 
obtained analytically from the strong coupling mapping 
of the ionic Hubbard model onto an effective Heisenberg 
model. 

We have also studied the bond-order parameter, the 
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order parameter associated with dimerization, as well as 
the associated bond-order susceptibility. The result of 
the delicate finite-size extrapolation indicates that there 
is a finite bond-order parameter in the intermediate re- 
gion between Ud and Uc2- There is a divergence in the 
bond-order susceptibility at both Ud and at Uc2, as one 
would expect from two continuous quantum phase tran- 
sitions. However, the bond-order susceptibility diverges 
in the entire strong coupling phase U > Uc2, albeit more 
weakly than at Ud- We have pointed out that this is in 
accordance with the behavior found in the strong cou- 
pling phase of the ordinary Hubbard model. 

We find that the electric susceptibility is finite for U < 
Ud but diverges roughly as L^-^ at Ud- This divergence 
is weaker than the one found for non-interacting electrons 
(with (5 = 0) and in the metallic phase of the t-i'-Hubbard 
modeli^ii A finite-size scaling analysis of both the bond- 
order susceptibility and the electric susceptibility yield 
the same critical exponents at Ud- However, the value, 
ryi ~ 0.45, is not consistent with the critical exponents of 
the classical two-dimensional Ising modeliS 

The electric susceptibility also seems to diverge, al- 
beit quite weakly, for U > Uc2- Correspondingly, the 
density-density correlation function has a long-distance 
decay which is of power-law form, but with a small pref- 
actor which becomes smaller with increasing J7, and a 
relatively large exponent of approximately 3 — 3.5. We 
speculate that this behavior is related to mixed spin and 



charge character of excitations present in the strong cou- 
pling phase of the ionic Hubbard model, in contrast to 
the ordinary Hubbard model. 

We point out that the divergence of the electric suscep- 
tibility at Ud and for U > Uc2 does not necessarily imply 
a finite Drude weight. Based on our results for various 
energy gaps and the electric susceptibility, we therefore 
cannot unambiguously classify all different phases and 
transition points as being metallic or insulating. 

Finally, we have presented DMRG results for the po- 
sition of the crossing of the ground state and the first 
excited J7x and the crossing of the first two excited states 
t/xx on systems with periodic boundary conditions on up 
to 64 sites. The finite-size extrapolation of J7x gives Ud- 
Due to the loss of accuracy, it is somewhat less clear that 
the finite-size extrapolation of ?7xx corresponds to Uc2 - 
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